format long;
clear;
L=5e-3;
C=50e-9;
R=1e3;
f=[0:100:500e3];
[h,v] = size(f);
hw = zeros(v);
hw(1,1) = 0;
hw(1,2) = 0;
for i=2:v
  CW = 1/(2*pi*f(i)*1j*C);
  LW = 2*pi*f(i)*1j*L;
  hw(i,1) = f(i);

  %r = (CW/(CW+LW+R));
  r = (LW/(CW+LW+R));
  hw(i,2) = abs(r);
  hw(i,3) = angle(r);
end;

fscan = [10000:200:50e3];
[h1,v1]=size(fscan);

for i=1:1
  fs = fscan(i) * 128;
  N = fs / (fscan(i) / 4);
  firN = N - 1;
  hnw(1, 1) = 0;
  hnw(1, 2) = 0;
  for j = 1: N
     nf = fscan(i) / 4 * j;
     hwfi = 0;
     hwti = 100;
     hwf = 1;
     hwt = 1;
     haf = 0;
     hat = 0;
     for k = 1:(v - 1)
        if((nf >= hw(k,1))&&(nf <= hw(k + 1, 1)))
          hwfi = hw(k, 1);
          hwf = hw(k, 2);
          haf = hw(k, 3);

          hwti = hw(k + 1, 1);
          hwt = hw(k + 1, 2);
          hat = hw(k + 1, 3);

          break;
        end
     end
     pk = ((hwt - hwf) / (hwti - hwfi));
     hwn = pk * nf+ hwf - pk * hwfi;


     if(hat < haf)
      hat = hat + 2 * pi;
     end

     pa = ((hat - haf) / (hwti - hwfi));
     han = pa * nf+ haf - pa * hwfi;

     hnw(j + 1, 1) = nf;
     hnw(j + 1, 2) = hwn;
     hnw(j + 1, 3) = han;
  end
  firn = calzfir(hnw(:,2), hnw(:,3), N, fs /N, 0);
  [h2,v2] = size(hnw);


  for m = 1: 1
    hn = m;
    ampt = 0;
    f0 = fscan(i) * hn;
    for k = 1:h2
      if(abs(f0 - hnw(k, 1)) < 1)
        ampt = hnw(k,2);
        break;
      end
    end
    ns = 1024;
    tfs = 1 / fs;
    pr = 1 / f0 / tfs;
    xs = 0:1:(ns - 1);
    xts = xs * tfs;


    sap = 1 / fscan(i) / tfs;
    for i=0:(ns - 1)
      ni = rem(i, sap);
      if(ni <= sap / 2)
         ydata(i + 1) = 1;
      else
         ydata(i + 1) = 0;
      end
    end
    %ydata = sin(2 * pi * f0 * xts);
    yfilter = zfilter(ydata, ns, firn, N);
    pr = pr * 4;
    yfn = yfilter(pr + 1 : ns);
    amp = max(yfn);
    fftyfn = fft(yfn);
    ni = (f0 / (fs / length(yfn)));

    nangle = angle(fftyfn(ni + 1));
    ampa(m, 1) = amp;
    ampa(m, 2) = nangle;
    ampa(m, 3) = ampt;
  end

  plot(ydata, "r");
  hold on;
  plot(yfn, "g");

  ylim([-1, 1]);
end

